I know I billed this as “Ditching ArcGIS” in the poll, but I 1) am not an ArcGIS expert and 2) there may still be things that you must do in ArcGIS to get your work done. If that’s the case, there is the R - ArcGIS Bridge that makes it possible for these two analysis tools to play with one another, BUT I won’t discuss that here…
# Set global knitr chunk options
if (.Platform$OS.type == "unix") {
# Do not specify Cairo device for MacOS
knitr::opts_chunk$set(warning = F, message = F)
} else {
knitr::opts_chunk$set(warning = F, message = F, dev.args = list(type = "cairo"))
}
You may need to install some packages (and their dependencies) that you don’t already have.
# List packages required to run the script -------------------------------------
pkgs <- c("tidyverse","ggmap","sf","cowplot","here",
"knitr","ggrepel","rnaturalearth","rgeos",
"FedData","raster")
# Install and load all CRAN packages provided from a character vector
load_pkgs = function(pkgs) {
new_pkgs = pkgs[!(pkgs %in% installed.packages()[ ,'Package'])]
if (length(new_pkgs) > 0) install.packages(new_pkgs,repos = "http://cran.cnr.berkeley.edu/")
invisible(lapply(pkgs,function(x)
suppressPackageStartupMessages(library(x,character.only = T))))
}
# Load packages
load_pkgs(pkgs)
dir.create(here("Figs"))
dir.create(here("Output"))
# Define lat and long bounds for west coast map
wc.lat <- c(32, 52)
wc.long <- c(-130, -116)
# Define lat and long bounds for Monterey Bay
mb.lat <- c(36.5, 37)
mb.long <- c(-122.2, -121.7)
# Read locations
locations <- read_csv(here("Data/locations.csv")) %>%
filter(group %in% c("city", "landmark") )
# Set west coast boundaries for stamen maps
wc.bounds.stamen <- c(left = min(wc.long), bottom = min(wc.lat),
right = max(wc.long), top = max(wc.lat))
# Download stamen map of west coast
wc.map.stamen.toner <- get_stamenmap(wc.bounds.stamen, zoom = 6, maptype = "toner-lite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw()
# Add layers to map
wc.stamen1 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text(data = locations, aes(long, lat, label = name)) +
ggtitle("Basic options")
wc.stamen2 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text(data = locations, aes(long, lat, label = name),
colour = "red", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) +
ggtitle("Formatted labels")
wc.stamen3 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
ggtitle("Formatted and repelled lables")
# Combine maps
wc.grid.stamen <- plot_grid(wc.stamen1, wc.stamen2, wc.stamen3, nrow = 1)
# Save map
ggsave(wc.grid.stamen, filename = here("Figs/wc_map_stamen_toner.png"), width = 12, height = 7)
# Print map
include_graphics(here("Figs/wc_map_stamen_toner.png"))
# Set west coast boundaries for stamen maps
wc.bounds.google <- c(mean(wc.long),mean(wc.lat))
# Download stamen map of west coast
wc.map.google <- get_map(location = wc.bounds.google, zoom = 4, source = "google", maptype = "satellite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
xlim(wc.long) + ylim(wc.lat)
wc.google1 <- wc.map.google +
geom_point(data = locations, aes(long, lat), colour = "white") +
geom_text(data = locations, aes(long, lat, label = name),
colour = "white", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) +
ggtitle("Minor formatting")
wc.google2 <- wc.map.google +
geom_point(data = locations, aes(long, lat), colour = "white") +
geom_text_repel(data = locations, aes(long, lat, label = name),
colour = "white", size = 2) +
ggtitle("Repelled labels")
# Arrange plots
wc.grid.google <- plot_grid(wc.google1, wc.google2, nrow = 1)
# Save map
ggsave(wc.grid.google, filename = here("Figs/wc_map_google.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_google.png"))
# Set west coast boundaries for stamen maps
mb.bounds.google <- c(mean(mb.long),mean(mb.lat))
# Download stamen map of west coast
mb.map.google <- get_map(location = mb.bounds.google, zoom = 10,
source = "google", maptype = "satellite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
xlim(mb.long) + ylim(mb.lat)
mb.google1 <- mb.map.google +
geom_point(data = filter(locations, name != "Monterey Bay"), aes(long, lat), colour = "white") +
geom_text(data = filter(locations, name != "Monterey Bay"), aes(long, lat, label = name),
colour = "white", size = 4, vjust = 0, nudge_y = 0.01)
# Save map
ggsave(filename = here("Figs/mb_map_google.png"), width = 7, height = 7)
# Print map
include_graphics(here("Figs/mb_map_google.png"))
rnaturalearth for shoreline and country files. Below is a comparison of the medium- and large-scale country datasets.
# Import NaturalEarth coastlines and countries
coast_sf <- ne_coastline(scale = "medium", returnclass = "sf")
coast_sf_lg <- ne_coastline(scale = "large", returnclass = "sf")
countries_sf <- ne_countries(scale = "medium", returnclass = "sf")
countries_sf_lg <- ne_countries(scale = "large", returnclass = "sf")
# View unique regions
# sort(unique(countries_sf$region_wb))
# sort(unique(countries_sf$subregion))
# Select only certain regions to speed plotting
na_sf <- filter(countries_sf, subregion %in% c("Northern America","Central America"))
na_sf_lg <- filter(countries_sf_lg, subregion %in% c("Northern America","Central America"))
# Create West Coast map (medium scale)
wc.ne <- ggplot() + geom_sf(data = na_sf) +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
scale_x_continuous(name = "Longitude", limits = wc.long) +
scale_y_continuous(name = "Latitude", limits = wc.lat) +
theme_bw() + coord_sf() +
ggtitle("Medium scale")
# Create West Coast map (large scale)
wc.ne.lg <- ggplot() + geom_sf(data = na_sf_lg) +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
scale_x_continuous(name = "Longitude", limits = wc.long) +
scale_y_continuous(name = "Latitude", limits = wc.lat) +
theme_bw() + coord_sf() +
ggtitle("Large scale")
# Create map grid
wc.grid.ne <- plot_grid(wc.ne, wc.ne.lg, nrow = 1)
# Save map
ggsave(wc.grid.ne, filename = here("Figs/wc_map_natEarth.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_natEarth.png"))
sf packageThe elevatr package provides easy access to elevation data from a variety of sources.
FedData Github repo and tutorial with examples of various data sources (e.g., .
Presently not executed; still under development
# Get a map of Moss Landing
ML.bbox <- get_map("Moss Landing, CA", zoom = 14)
# Plot the map
ggmap(ML.bbox)
# Get the bounding box from its attributes
bb <- attr(ML.bbox, "bb")
# Create an extent polygon from the Monterey Bay bounding box
ML.extent <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon,
bb$ll.lat, bb$ur.lat),
proj4string="+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
# Map the bounding box to check the results
ggmap(ML.bbox) + geom_polygon(aes(x=long, y=lat), size=3, color='purple',
data=ML.extent, fill=NA)
# Download National Elevation Data
ned_ML <- get_ned(template = ML.extent, label="ned_ML",
res="1", force.redo = F, raw.dir = here("Data"))
# Convert raster to data frame for plotting in ggplot
ned_df <- rasterToPoints(ned_ML) %>%
data.frame() %>%
rename(lat = y,
long = x,
elev = grdn37w122_1)
# Map elevation
ggplot(ned_df, aes(long,lat,fill = elev)) + geom_raster() +
geom_contour(aes(z = elev), alpha = 0.5, binwidth = 1) +
scale_fill_viridis_c(name = "Elevation (m)") +
theme_bw()
Simple Feaures for R: Blogs, presentations, vignettes for the sf package
A Tidy Approach to Spatial Data: Simple Features from our own Eric Anderson
GIS with R: an introduction by Francisco Rodriguez-Sanchez (Pakillo)
An Intro to GIS with R by Jess Sadler (covers sp, sf, and rnaturalearth, among other things).